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We investigate the effect of eiectron-phonon interaction on the phononic properties in the one- 
dimensional half-filled Holstein model of spinless fermions. By means of determinantal Quantum 
Monte Carlo simulation we show that the behavior of the phonon dynamics gives a clear signal 
of the transition to a charge-ordered phase, and the phase diagram obtained in this way is in 
excellent agreement with previous DMRG results. By analyzing the phonon propagator we extract 
the renormalized phonon frequency, and study how it first softens as the transition is approached 
and then subsequently hardens in the charge- ordered phase. We then show how anharmonic features 
develop in the phonon propagator, and how the interaction induces a sizable dispersion of the dressed 
phonon in the non-adiabatic regime. 



INTRODUCTION 

Despite enormous theoretical efforts, effects arising 
from the eiectron-phonon interaction in strongly corre- 
lated many-body systems still remain incompletely un- 
derstood. Achieving an understanding of the complicated 
interplay between electrons and phonons is, however, 
essential to explain such diverse phenomena as colos- 
sal magnetoresistance in the manganitespj, the Peierls 
instability in quasi-lD materials [2j, and high tempera- 
ture superconductivity in alkali- metal doped fullerenesQ 
and cuprate compounds 4] . Even in moderately or 
weakly correlated materials, the eiectron-phonon inter- 
action can give rise to interesting effects which are not 
understandable in the framework of the standard theo- 
ries of eiectron-phonon interaction, namely the Migdal- 
Eliashberg theory of superconductivity and the Born- 
Oppenheimer adiabatic principle. As a notable example, 
we mention the anomalies in the phononic properties of 
superconducting MgB2 0- 

In this work we present a numerical investigation of 
a simple model for coupled eiectron-phonon systems — 
the half-filled Holstein model for spinless fermions — and 
mainly focus on the dynamics of the phononic degrees of 
freedom. In spite of its apparent simplicity, the model 
fully accounts for the competition between local quantum 
fluctuations and the tendency for charge ordering, and is 
thus a powerful tool to understand the physics of more 
realistic systems. In one-dimension the system forms a 
metallic Luttinger liquid for weak eiectron-phonon cou- 
pling. As the coupling is increased, a quantum phase 
transition occurs to an insulating state with long-range 
charge-density-wave (CDW) order. A wide variety of 
methods have been used to investigate this phase tran- 
sition, with varying degrees of success. Exact diagonal- 
ization schemes |6j permit extremely accurate calcula- 
tion of the ground-state and low-lying excitations, but 
are limited to treating rather small clusters due to the 
large Hilbert space required for the phonon degrees of 



freedom, and are thus subject to large finite-size effects. 
Larger systems are accessible using the density matrix 
renormalization group (DMRG) method, and in a recent 
work p| Bursill et al were able to determine the location 
of this phase boundary using this approach. Their result 
correctly recovers the adiabatic and anti-adiabatic limits 
(see FigflJ which can be evaluated analytically ||, and 
interpolates smoothly between them. Results obtained 
from a variational Lanczos scheme on small clusters @ 
agree well with the DMRG result, but data produced 
by Quantum Monte Carlo (QMC) methods, such as the 
worldline QMC method used in the pioneering investi- 
gation of Ref.|8| and Green's function QMC 10], show 
significant deviations, the cause of which is not known. 
The role of quantum lattice fluctuations on observables 
has been discussed in Refs. 0,0], while the spectral 
properties have been analysed in Refs. (l^fl4| 

In this work we make use of the determinantal QMC 
method |15| which allows us to compute dynamical prop- 
erties like the electron and phonon propagators without 
approximations. These quantities are not equally acces- 
sible to DMRG, which is otherwise the best method to 
establish the phase diagram of one-dimensional models. 
In contrast to previous QMC investigations we are able 
to obtain estimates for the location for the phase bound- 
ary in excellent agreement with the "benchmark" DMRG 
result. We concentrate on the phonon propagator, which 
appears to produce more robust results than the elec- 
tronic properties considered in previous QMC studies, 
and show firstly how its behavior reveals the location 
of the phase-boundary. We then go on to investigate 
in detail the phonon-softening and anharmonicity effects 
induced by the eiectron-phonon interaction by combin- 
ing analysis of the imaginary-time data with the spectral 
functions obtained by analytically continuing this data 
to real frequencies. This allows us to identify a renor- 
malized phonon frequency, and to follow its behavior as 
a function of eiectron-phonon coupling. Beside a generic 
softening of the frequency, we observe how for large bare 



2 







1 | 1 III 








□ WL-QMC 








a GF-QMC 


7 






« <> DMRG 


/ 

7 


3 




. D-QMC 

mean field 

■ — strong coupling 


// - 

A' 

// 


3° 

- 2 
"so 


' v.. 

A • ^ 


LDW insulator ,-■ ' 

o' / 






• N 

\ ° 
*\ 


» V 

□ 


1 






V 







Luttinger liquid metal • . 

1 1 i i 




0.1 


1 10 



t/co 



FIG. 1: Phase diagram of the ID Holstein model, showing 
the phase boundary between the CDW insulator and Lut- 
tinger liquid metal. For ease of comparison with previous re- 
sults, the y-axis is in units of g' = g/\/2mui, the coupling be- 
tween the electron and the quantized phonon operators a% /a\ . 
Squares indicate results from the world-line QMC method 
from Ref.Q, triangles are results from the Green's function 
QMC investigation of Ref. llCll and the solid circles are the 
results of this work using determinantal QMC. Diamonds de- 
note the results of the DMRG investigation in Ref.0. 

phonon frequency, the harmonic approximation for the 
dressed phonon propagator fails as soon as the coupling 
with the electrons is introduced. We also discuss how the 
renormalization of the phonon properties also gives rise 
to a phonon dispersion — even if the bare phonons are 
dispersionless. 

The paper is organized as follows. In Sec. II we intro- 
duce the model and some details about the QMC sim- 
ulations. Sec. Ill presents the results, and it is in turn 
divided into two subsections — the phonon propagator 
in imaginary time and its analytical continuation to real 
frequencies. In Sec. IV we summarize and give our con- 
clusions. 



Here c i /c i are the fermion annihilation/creation opera- 
tors and n, is the fermion number operator cjc^. The 
electron-phonon coupling is set by the phonon fre- 
quency is given by ojq, and <ji and pi denote the phonon 
displacement and momentum operators respectively. We 
will express all energies in units of the electronic hopping 
t, and set the phonon mass m equal to one. Unlike Ref.fg 
we work in the grand-canonical ensemble, in which the 
electronic density is regulated by the chemical potential 
fj,. In this work we only consider the case of the half-filled 
system, which is given by /i = g 2 /2uJq. 

To simulate this model we employ the well-known 
determinantal QMC method (DQMC) developed by 
Blankenbecler, Scalapino and Sugar 15]. In this ap- 
proach the fermion degrees of freedom are analytically 
integrated out of the action, which is straightforward 
as the Hamiltonian Q is bilinear in fermion operators, 
leaving an effective action expressed just in terms of the 
phonon displacement field. By formally replacing the 
time coordinate with imaginary-time (r = it), the par- 
tition function of the model can then be simulated as a 
path integral of a Euclidean field theory using standard 
Monte Carlo techniques. In this formalism the T-axis 
represents an additional compact dimension, the extent 
of which is given by the inverse temperature 0. In or- 
der to sample the zero-temperature properties of the sys- 
tem, it is important that sufficiently large values of 
are used. By comparing the convergence of simulations 
as was increased we established that in the adiabatic 
regime (u>q < t) it was sufficient to take 0t = 8, but for 
the highest phonon frequency that we studied (u>o = it) a 
larger value of 0t = 16 was required. Increasing the size 
of also requires increasing the number of time-slices 
used in the simulation, in order to keep the systematic 
error arising from the Trotter decomposition sufficiently 
small. This both increases the simulation's running time 
and also diminishes its numerical stability, and therefore 
sets an upper limit on the value of coq we are able to treat 
with this method. 

Observables, such as the phonon correlation function, 
are obtained from the simulation as thermal averages of 
the form: 



MODEL AND SIMULATION 

We consider the following Hamiltonian, describing 
spinless fermions moving on a periodic ID lattice and 
interacting with a dispersionless phonon at each lattice 
site: 

H = - t^T (c\c l+1 + ff.c.J -9^2 n iH +fJ-^2n % 



Ai(r) = ( qi (T) qj (0)), O<r<0 

= Tr[ qi (r) qj (0)e^ H ]/Z (2) 

where Z = Tr[exp(— 0H)], and time-dependent opera- 
tors arc defined as q(r) = cxp(— Hr) q exp(HT). It is im- 
portant to note that the correlation functions produced 
by the simulation depend on the imaginary, or Matsub- 
ara, time-coordinate r. To study the system's dynamical 
properties it is thus necessary to make an analytic contin- 
uation of these functions to the real-time domain. This 
amounts to the solution of the following inverse problem 



3 



0: 



S(r) 



duj- 



-0W 



XtM, 



(3) 



where S(t) is the Matsubara correlation function, and 
Xt( w ) ^ s the imaginary component of the time-ordered 
susceptibility. To perform this inversion we use a max- 
imum entropy technique (lrj based on the singular- 
value decomposition |l8| of the kernel of Eq|3J By im- 
posing positivity and smoothness constraints on the so- 
lution we find that the instabilities typically associated 
with numerical analytical continuation procedures can be 
controlled, and this method is thus able to produce stable 
results of high resolution for the phonon spectral func- 
tion. 



RESULTS 

The Matsubara phonon propagator 

In principle, the transition from the Luttinger liquid 
phase to the CDW insulator can be observed in the 
staggered-phonon order parameter 
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(4) 



which takes a non-zero value in the CDW phase. It is 
problematic, however, to measure this quantity directly 
in a simulation of a finite system, as in the CDW phase 
the system continually tunnels between the two degen- 
erate ground-states, causing the expectation value of m p 
to vanish in the ergodic limit. Over short runs, however, 
a non-zero value can be obtained for m p if the system 
remains trapped in one of the minima for the duration 
of the simulation. This allows the location of the phase 
transition to be located approximately but this esti- 
mate is intrinsically rough, and cannot be improved easily 
as the effect of increasing the number of measurements 
is to make m p vanish. Nonetheless it provides a useful 
initial comparison for other estimates, and we found that 
in all cases it was consistent with the values we obtained 
by other methods. 

As mentioned above, our interest is devoted to the 
phonon dynamics, which provides us with an alternative 
way to pinpoint the phase transition. More precisely, 
we study the phonon propagator @, which in the non- 
interacting case (g = 0) reads: 



D?,-(t,Wo) = 



1 cosho;o(r-/3/2) 
2luq sinhwo/3/2 



(5) 



It should be noted that as the Holstein phonon is purely 
local, the bare phonon propagator also only has a trivial 
spatial dependence. For small electron-phonon coupling 



it is reasonable to expect that the phonon propagator has 
the same form, but with with a renormalized (softened) 
frequency SI given by S1 2 /u)q = 1 — U(0)/lu 2 , where II(cj) 
is the local phonon self-energy in real frequencies. This 
ansatz is most accurate when the frequency-dependence 
of the phonon self-energy is weak, and deviations from 
this ideal behavior are the fingerprints of phonon an- 
harmonicity, i.e. of the difficulty in describing the fully 
dressed phonon as a single renormalized oscillator. 

As the system approaches the CDW transition, the 
staggered (k = tt) phonon correlation function is ex- 
pected to soften and eventually produce a sharp peak at 
zero frequency, the weight of which measures the conden- 
sation of the phonons. For an infinite system at T = 0, 
this quantity would be the phonon staggered-order pa- 
rameter m p . On the other hand, the appearance of zero- 
frequency weight may also occur in the local (k = 0) 
correlation function, which is not, however, directly as- 
sociated with the CDW transition. In this case, the shift 
of weight can be associated with a polaron crossover from 
a good metal to a bad metal, in which the electronic mo- 
bility is strongly reduced by the large coupling to the 
phonons. 

We tested this expectation by fitting the local and stag- 
gered imaginary time correlation functions with the func- 
tion D°(t, fl) + c, using the softened frequency and the 
shift c as fitting parameters. The rigid shift in imagi- 
nary time c is associated with a 5-like peak at zero fre- 
quency, and describes the static average of the phonon 
field. For the local correlator the value of c indicates the 
static uniform distortion, which reduces the electron mo- 
bility through polaronic effects, while for the staggered 
case it is simply equal to the staggered phonon order 
parameter m p . For a small value of the bare phonon- 
frequency ujq = 0.5, we find that using the above form 
for the fitting function yields excellent results for both 
the local and staggered propagators, at all values of the 
electron-phonon coupling. We show one such example for 
the staggered propagator in FigJ^,. The good quality of 
the fit is corroborated by the form of the Fourier trans- 
form of the correlator, shown to the right. The Fourier 
representation of the non-interacting propagator (jSJ) can 
be easily shown to be 



D°Jiu; n ,ujo) 



(6) 



where u n are the Matsubara frequencies, Lu n = 2im//3. 
Simply replacing the bare frequency, ujq, in Eq0 with 
the renormalized frequency, f2, obtained from the fit- 
ting procedure provides an extremely good fit to the 
Fourier transform of the Matsubara data, as can be seen 
in Fig[2l3. A similar behavior is observed for u>o = 1 
(Figs|2t,d), where the fitting procedure again yields ex- 
cellent results for both the Matsubara-time propagator 
and its Fourier transform. We emphasize that these fits, 
which essentially pass through every data-point, are ob- 



FIG. 2: (a) The staggered phonon propagator as a function 
of imaginary time for luq — 0.5, g — 0.6 (c) for u>o — 1 and 
g = 1.5 and (e) for wo = 4 and g = 14.0. All data is taken from 
simulations of 16-site systems. The line shows the fit to the 
form D(t, u>o) = D°(t, Q)+rrip. In (b), (d) and (f) we compare 
the corresponding Fourier transforms of these propagators in 
Matsubara frequency space with the renormalized Lorentzian 
Eq© 

tained using merely two fitting parameters with a direct 
and suggestive physical meaning. 



On increasing the bare phonon frequency further, how- 
ever, we find that the fitting function no longer re- 
produces the data once the electron-phonon coupling is 
turned on. As can be seen from Fig|2? the curvature 
of the phonon propagator as a function of r cannot be 
described simply in terms of a softened frequency. This 
effect can be seen equally distinctly in frequency space 
in Fig|2k, where D(iuj n ) clearly deviates from the pre- 
dicted Lorentzian behavior. Effectively the non-adiabatic 
bare phonon introduces anharmonic effects, which are en- 
hanced by increasing the coupling. The relationship be- 
tween non-adiabaticity and anharmonicity has been pre- 
viously proposed for MgB 2 ll9|. 

In FigOHwe show the values for the renormalized fre- 
quencies, obtained by the fitting procedure, for both the 
local and staggered phonon propagators. The values of 
rrip produced by the fitting method are also shown below. 
For ujq = 0.5 (Figjjjji) it can be seen that in the metal- 
lic region the local phonon frequencies smoothly reduce 
from their bare values as g increases, to reach a minimum 
value at the point g/u>o = 1.6. Throughout this process 
the staggered phonon-frequency is softened to a greater 
degree than the local phonon. At the local minimum, m p 
suddenly acquires a non-zero value, indicating that this 
point signals a transition to the CDW regime, and thus 
most of the weight in the phonon propagator moves to 
the zero frequency peak. From Fig^ it can be clearly 
seen that this estimate compares extremely well with the 



FIG. 3: Softened phonon frequency, Q, and order parameter, 
mp, obtained by fitting the phonon propagator to the form 
D(r) = D°(t, $7) + nip. Solid lines indicate the local phonon 
frequency, dotted lines the staggered phonon. (a) uiq = 0.5, 
(b) ujo = 1 and (c) ujq — 4. In all cases the staggered phonon 
is softened more than the local phonon, and it can be clearly 
seen that the minimum becomes increasingly deep (i.e. the 
phonon softens more) as ujq is increased. 

location of the phase transition found in Ref.Q- On in- 
creasingo further, the renormalized frequency is seen to 
harden |20( , and eventually approach the bare frequency 
in the atomic limit, where the system is no longer metal- 
lic and therefore there is no screening of the phonons. 

For higher phonon-frequencies a similar behavior 
occurs, with the softened phonon frequency passing 
through a well-defined minimum at which turns on. 
For ujq — 1 (FiglJb) the softening of the phonon fre- 
quency is more pronounced, with the staggered phonon 
frequency decreasing by a factor of one-half as compared 
with a factor of about two-thirds for the previous case. 
This trend continues for the antiadiabatic case, ujq = 4 
(FigEt) , for which the bare phonon frequency is reduced 
by almost a factor of ten. Since the fitting procedure 
is less trustworthy for this case, the values obtained 
for the fitting parameters should be treated more cau- 
tiously. The good agreement with the DMRG result that 
is nonetheless obtained indicates that the procedure is 
successfully describing the phonon dynamics with rea- 
sonable accuracy, and underlines the reliability of using 
this behavior as a signal of the phase transition. 



Further insight about the effect of the interaction on 
phononic properties can be gained by the analysis of the 
statistical distribution of the phonon displacement field 
qj. In Fig. 0]we show how this quantity evolves as the 
electron-phonon coupling is increased. A clear qualita- 
tive change occurs from an unimodal distribution in the 
metallic phase, in which q — is the most probable value, 
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FIG. 4: Histograms of the phonon displacement field in 
the vicinity of the phase transition for a system with phonon 
frequency Wq = 1. (a) g = 0.0, (b) g = 2.25 and (c) g = 2.5. 
Note how the non-interacting Gaussian distribution centered 
on zero first flattens as the phase transition is approached, 
and then splits into two in the CDW regime, indicating the 
formation of a polaronic state. 

to a bimodal distribution with maxima |2l| separated by 
2<Zo = <?/wq, which signals the formation of local lattice 
distortions peaked around ±qq. In previous work [2^ . 
the evolution of the distribution of the displacements 
has been used to characterize the polaron crossover in 
the Holstein model. In principle, the formation of local 
distortions does not automatically imply a CDW order- 
ing, as has been discussed for example in Ref.j2^]. For 
the system we study, however, it does appear that the 
formation of local lattice distortions and CDW ordering 
do occur simultaneously, implying that a polaronic metal 
state is not present in the transition region. 



Analytic continuation of the phonon propagator 

To complement the previous analysis of the Matsubara 
correlation functions, we now study the real-frequency 



propagator by employing a maximum entropy method 
to make the analytic continuation from Matsubara time 
to real-frequencies. This avoids the need of assuming a 
given analytic form for the propagator, and also allows 
us to study the softening of the phonon in more detail. 
In Fig|S] we present contour plots of the local phonon 
spectral function for different bare phonon-frequencies. 
Darker areas correspond to larger weights. In all figures, 
the weight is clearly concentrated at the bare frequency 
(±Cl>o) at g = 0, and for weak coupling it remains concen- 
trated in a single feature at the renormalized frequency 
il. As the renormalized phonon softens further, a zero 
frequency peak appears and takes most of the weight. 
In terms of the fitting parameters, the appearance of this 
feature corresponds to a non-zero value of the shift c, im- 
plying a sizable static lattice deformation, which reduces 
the mobility of the electrons. After the CDW transition 
a higher energy phonon branch forms, and, as was seen 
previously from the Matsubara analysis, its energy hard- 
ens and eventually converges to the bare frequency in the 
limit of extreme strong-coupling. In Fig|S]the formation 
of this branch is clearly visible for the non-adiabatic cases 
u>o = 1 and 4, and although it also occurs in the adia- 
batic case the higher energy features are rather obscured 
by the large zero frequency peak. Thus, in the limit of 
strong-coupling the system recovers the appearance of 
the atomic limit in which the phonon is unrenormalized, 
due to the lack of metallic screening. 



A similar analysis can be carried out for all the val- 
ues of the exchanged phonon momentum k, resulting in 
a momentum dependent renormalized phonon frequency 
fi(fc). We plot this quantity for u>o = 1 and g = 1 and 
2 in Fig. [BJ In the weaker coupling case, the renor- 
malized phonon is still basically dispersionless like the 
bare Holstein phonon. For the larger coupling, however, 
the renormalized phonon becomes dispersive, exhibiting 
an approximately cos(fc/2) dependence on momentum. 
This dramatically demonstrates that if a material has a 
sizable electron-phonon interaction, the bare dispersion 
of the phonons may be substantially different from the 
fully dressed one. This has to be seriously taken into ac- 
count in the derivation of effective models, in which the 
experimentally observed fully dressed dispersion should 
not be used as a bare dispersion for a model calculation 
to avoid double counting. 



CONCLUSIONS 



In this paper we have investigated the effect of 
electron-phonon interaction on the half-filled one- 
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FIG. 5: Analytic continuation of the local phonon propagator Du(t) for phonon frequencies luq = 0.5, 1.0,4.0. 
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FIG. 6: Phonon susceptibility \" (k,u)) for coo — 1. The peaks 
in this quantity (dark areas) indicate the dispersion relation 
followed by the dressed phonon. Left: g = 1.0, the dressed 
phonon is practically dispersionless, resembling the bare Hol- 
stein phonon. Right: g = 2.0, the phonon peak shifts toward 
zero as k increases indicating that the dressed phonon has 
acquired a substantial dispersion from the electron-phonon 
interaction. 



dimensional Holstcin model by means of DQMC simu- 
lations. We have firstly shown how examination of the 
phonon propagator can be used to determine the bound- 
ary for the CDW phase transition, which agree well with 
the state-of-the-art DMRG data. 

The agreement of our results with the DMRG phase 
diagram makes us confident about the reliability of our 
method in the evaluation of other observables, which are 
more difficult to obtain with DMRG, such as dynamical 
properties. 

In particular, we focus on the phonon propagator, in 



order to discuss the dressing of the phonon degrees of 
freedom for large coupling with the electrons. An analy- 
sis of the Matsubara phonon propagator allows a renor- 
malized phonon frequency to be deduced directly from 
the data produced by the QMC simulation. This has 
revealed how the phonon frequency softens as the cou- 
pling is increased until, at a critical value of <?, the order 
parameter m? v turns on and the system makes a phase 
transition to the CDW regime. The degree of softening is 
considerably higher in the non-adiabatic regime as com- 
pared to the adiabatic case. Our analysis has also high- 
lighted the relation between non-adiabatic effects and an- 
harmonicity. When the bare phonon frequency is small, 
the phonon propagator describes basically a single har- 
monic phonon even for large coupling, while in the non- 
adiabatic case, the electron-phonon coupling induces an- 
harmonic effects through a frequency dependent phonon 
self-energy. 

By using analytic continuation methods we were able 
to study the phonon renormalization in detail. At the 
critical point we have seen how the phonon mode splits 
into two — one mode developing into a soft mode, coex- 
isting with the other mode which then hardens and re- 
approaches the bare frequency as g is increased further. 
This technique has also revealed the strong momentum- 
dependence of the dressed phonon in the non-adiabatic 
regime, arising from the electron-phonon interaction. In 
particular, a sizable phonon dispersion arises as the cou- 
pling is increased, even starting from a bare dispersionless 
phonon, suggesting that care must be taken in building 
up models for electron-phonon interaction starting from 
experimentally observed phonon properties. 
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